% Norm

function y = func_norm(m, CONSTS)

    c = CONSTS.c;
    a = CONSTS.a;
    
    r_high = 1000;
    num_r_points = 1000;
    r_step = r_high/num_r_points;
    r1 = (0 : r_step : 1)';
    r2 = (1 : r_step : r_high)';
    
    inside  = func_inside(r1, m, CONSTS);
    outside = func_outside(r2, m, CONSTS);
    
    inside(isnan(inside)) = 0;
    outside(isnan(outside)) = 0;
    
%     y = c*(a^2) * (trapz(r1, inside) + trapz(r2, outside));
    y = c * (a^2) * (quadgk(@(r)func_inside(r, m, CONSTS), 0.0, 1) + ...
            quadgk(@(r)func_outside(r, m, CONSTS), 1, 1e+20));
    
end

function y = func_inside(r, m, CONSTS)

    [E_rho_inside,   E_phi_inside,   foo]  = func_E_inside (r, m, CONSTS);
    [H_rho_inside,   H_phi_inside,   foo]  = func_H_inside (r, m, CONSTS);

    y = (E_rho_inside .* H_phi_inside + E_phi_inside .* H_rho_inside) .* r;

end

function y = func_outside(r, m, CONSTS)

    [E_rho_outside,  E_phi_outside,  foo] = func_E_outside(r, m, CONSTS);
    [H_rho_outside,  H_phi_outside,  foo] = func_H_outside(r, m, CONSTS);
    
    y = (E_rho_outside .* H_phi_outside + E_phi_outside .* H_rho_outside) .* r;

end

